About Issues Tutorials Documentation </span>

Example 1: Build and simulate DNA

This notebook builds a small DNA double helix, assigns a forcefield to it, and runs a molecular dynamics simulation.

Author: Aaron Virshup, Autodesk Research
Created on: July 1, 2016
Tags: DNA, molecular dynamics



In [ ]:

    
import moldesign as mdt
from moldesign import units as u

%matplotlib inline
from matplotlib.pyplot import *

# seaborn is optional -- it makes plots nicer
try: import seaborn  
except ImportError: pass

1. Create a DNA helix
2. Forcefield
3. Constraints
4. MD Setup
5. Minimization
6. Dynamics
7. Analysis

1. Create a DNA helix



In [ ]:

    
dna_structure = mdt.build_dna_helix('ACTGACTG', helix_type='b')
dna_structure.draw()



In [ ]:

    
dna_structure

2. Forcefield

The cell below adds forcefield parameters to the molecule.

NOTE: This molecule because is not missing expected atoms. If your molecule is missing atoms (e.g., it's missing its hydrogens), use the Forcefield.create_prepped_molecule method instead of Forcefield.assign.

Click on the ERRORS/WARNING tab to see any warnings raised during assignment.



In [ ]:

    
ff = mdt.forcefields.DefaultAmber()
ff.assign(dna_structure)

3. Constraints

This section uses an interactive selection to constrain parts of the DNA.

After executing the following cells, click on the 3' and 5' bases:



In [ ]:

    
rs = mdt.widgets.ResidueSelector(dna_structure)
rs



In [ ]:

    
if len(rs.selected_residues) == 0:
    raise ValueError("You didn't click on anything!")
    
rs.selected_residues



In [ ]:

    
for residue in rs.selected_residues:
    print('Constraining position for residue %s' % residue)
    
    for atom in residue.atoms:
        dna_structure.constrain_atom(atom)

Of course, fixing the positions of the terminal base pairs is a fairly extreme step. For extra credit, see if you can find a less heavy-handed keep the terminal base pairs bonded. (Try using tab-completion to see what other constraint methods are available)

4. MD Setup

This section adds an OpenMM energy model and a Langevin integrator to the DNA.



In [ ]:

    
dna_structure.set_energy_model(mdt.models.OpenMMPotential,
                               implicit_solvent='obc')

dna_structure.set_integrator(mdt.integrators.OpenMMLangevin,
                             timestep=2.0*u.fs,
                             temperature=300.0*u.kelvin,
                             frame_interval=1.0*u.ps)

You can interactively configure these methods:



In [ ]:

    
dna_structure.configure_methods()

5. Minimization

Nearly every MD simulation should be preceded by an energy minimization, especially for crystal structure data. This will remove any energetically catastrophic clashes between atoms and prevent our simulation from blowing up.



In [ ]:

    
trajectory = dna_structure.minimize(nsteps=200)
trajectory.draw()



In [ ]:

    
plot(trajectory.potential_energy)

xlabel('steps');ylabel('energy / %s' % trajectory.unit_system.energy)
title('Energy relaxation'); grid('on')

6. Dynamics

We're ready to run 25 picoseconds of dynamics at room temperature (that's 300º Kelvin). This will probably take a few minutes - if you're on an especially pokey computer, you might want to reduce the length of the simulation.



In [ ]:

    
traj = dna_structure.run(run_for=25.0*u.ps)
traj.draw()

7. Analysis

The trajectory object (named traj) gives direct access to the timeseries data:



In [ ]:

    
plot(traj.time, traj.kinetic_energy, label='kinetic energy')
plot(traj.time, traj.potential_energy - traj.potential_energy[0], label='potential_energy')
xlabel('time / {time.units}'.format(time=traj.time))
ylabel('energy / {energy.units}'.format(energy=traj.kinetic_energy))
title('Energy vs. time'); legend(); grid('on')



In [ ]:

    
# Using the trajectory's 'plot' method will autogenerate axes labels with the appropriate units
traj.plot('time','kinetic_temperature')
title('Temperature'); grid('on')

This cell sets up an widget that plots the RMSDs of any selected group of atoms. Select a group of atoms, then click "Run plot_rmsd" to generate a plot



In [ ]:

    
from ipywidgets import interact_manual
from IPython.display import display

rs = mdt.widgets.ResidueSelector(dna_structure)
def plot_rmsd(): 
    plot(traj.time, traj.rmsd(rs.selected_atoms))
    xlabel('time / fs'); ylabel(u'RMSD / Å')
interact_manual(plot_rmsd, description='plot rmsd')
rs



In [ ]: